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ABSTRACT 

We present in this paper both a linear study and numerical relativistic MHD simula- 
tions of the non-resonant streaming instability occurring in the precursor of relativistic 
shocks. In the shock front restframe, we perform a linear analysis of this instability in 
a likely configuration for ultra-relativistic shock precursors. This considers magneto- 
acoustic waves having a wave vector perpendicular to the shock front and the large 
scale magnetic field. Our linear analysis is achieved without any assumption on the 
shock velocity and is thus valid for all velocity regimes. In order to check our cal- 
culation, we also perform relativistic MHD simulations describing the propagation 
of the aforementioned magneto-acoustic waves through the shock precursor. The nu- 
merical calculations confirm our linear analysis, which predicts that the growth rate 
of the instability is maximal for ultra-relativistic shocks and exhibits a wavenumber 

1/2 

dependence oc k x . Our numerical simulations also depict the saturation regime of 
the instability where we show that the magnetic amplification is moderate but never- 
theless significant (8B/B ^ 10). This latter fact may explain the presence of strong 
turbulence in the vicinity of relativistic magnetized shocks. Our numerical approach 
also introduces a convenient means to handle isothermal (ultra-)relativistic MHD con- 
ditions. 

Key words: Shock waves - Magnetohydrodynamics (MHD) - Cosmic rays - Rela- 
tivistic processes 



1 INTRODUCTION 

Astrophysical shock waves and particle acceleration up to 
the highest energies are intimately connected processes. In 
non-relativistic supernova remnants a large amplification of 
the magnetic field at the forward shock has been detected 
by recent high angular resolution observations from X-ray 



satellites Chandra and XMM Newton (Ballet 20061. The 



very nature of the amplification process has been the object 
of intense research work in the last decade (see a recent 



review by Schure et al ( 2012 1) . It appeared that the current 



generated by the energetic particles (hereafter abusively 
considered as cosmic rays) streaming ahead of the shock 
front in the so-called CR precursor could provide such 
an amplification (IBelll 120041) even if the definite value of 
the amplified magnetic field remains uncertain ( |Riquelme| 
& Spitkovsky 2009[ ). In relativistic shocks, notably in 



gamma-ray bursts we can only indirectly evaluate magnetic 
field amplification in the shock precursor. Here again 
high magnetic amplification seems to be required in the 
upstream medium (Li & Waxman 2006). The fraction 



e_e of the shock energy flux imparted into the magnetic 
field has been inferred from multi-wavelength observations 
of gamma-ray burst afterglows and can reach values up 
to 1% (however see the discussion in |Lemoine| ( |2013[ ) ) . 
Here again the streaming of cosmic rays (CRs) ahead 
of the shock may provide the main source of free energy 
to generate the magnetic field ( |Milosavljevic fc Nakar|2006| ). 



In a previous work Pelletier et al ( 2009 1 (hereafter 
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PLM09) have investigated in detail the amplification of 
magnetic disturbances at relativistic shock waves. This 
results in two necessary conditions to get the Fermi 
acceleration process working in such environments. The 
first condition states that the turbulence triggered around 
the shock must be at small scale (with respect to the 
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triggering particle Larmor radius). The second condition 
states that the ratio of the perturbed to background 
magnetic field SB /Bo 3> 1 consistent with strong magnetic 
field amplification. PLM09 did perform 1 and 2D analytical 
analysis of magnetohydrodynamic (MHD) waves evolution 
in the CR precursor of highly relativistic perpendicular 
shocks. They have shown that Alfven waves are stable and 
that magneto-acoustic waves can be destabilized by the 
charge of CR ahead the shock. Hence it has been further 
shown that relativistic shocks can trigger microscopic 
instabilities that lead to turbulence at scales beyond the 



MHD regime (Lemoine fc Pelletier 2010 Bret 20091, one 



may also see the recent review by Bykov et al (20121. 
However, MHD waves are of great interest because they 
can provide a better confinement of high energy particles. 
They also can at least partly explain values of e_e inferred 
from observations. Finally they can also generate shock 
trains upstream the shock front that can contribute to 
particle acceleration and plasma heating. In this work we 
continue the analysis initiated by PLM09 by testing these 
analytical results with respect to relativistic MHD (RMHD) 
numerical simulations and by providing a more general 
linear analytical framework that can be generalized to more 
complex flow configurations. In this work we also provide 
an analysis of mildly- and non-relativistic flow regimes. 

The paper is organized as follows. In section [2] we dis- 
cuss the physical case considered in this work as well as the 
formalism used in the linear instability analysis. In section 
[3] we give a short presentation of the relativistic version of 
MPI-AMRVAC and hence describe the different stages of 
the streaming instability in a reference case corresponding 
to an ultra-relativistic shock wave. We end this section with 
a parametric survey of the instability with respect to the 
shock Lorentz factor and to the perturbation wavenumber. 
We conclude and provide a short discussion in section [4] 



2 NON-RESONANT STREAMING 

INSTABILITY NEAR RELATIVISTIC 
SHOCKS 

As pointed out by PLM09, ultra-relativistic astrophysical 
shock waves are likely to exhibit a perpendicular large scale 
magnetic structure, which means a magnetic field parallel to 
the shock front (hence perpendicular to the shock normal). 
Indeed, the Lorentz transformation of a magnetic field from 
the observer rest-frame (denoted with "obs" subscripts) to 
the shock rest-frame (denoted as "sh" subscripts) leads to 



B 
B, 



||, ah 



-B|| j0 b s 

r sh (Bj 



- A* x Ej 



(1) 



where /3 s h is the velocity of the shock wave in the observer 
frame (r^ 2 = 1 — /3 2 h ) while Ei )0 b s is the electromotive 
field occurring in the observer frame. Ultra-relativistic shock 
waves do have generically a perpendicular magnetic field in 
the frame of the shock, except if the magnetic field in the 
observer frame is aligned with the normal direction to the 
shock front within 0(l/r s h). This work considers hereafter 
the case of a pure perpendicular magnetic field in the shock 
frame where the magnetic field is oriented in a y direction 
parallel to the shock front. The shock front will be assumed 



as planar thus leading to fluid properties varying only along 
the normal to the shock front in the x direction. 



2.1 Special relativistic magnetohydrodynamics 
equations 

In this work we present time dependent numerical simula- 
tions of the evolution of MHD waves that propagate through 
the precursor of a perpendicular magnetized shock wave. 
The special relativistic magnetohydrodynamic (SRMHD) 
equations ruling both the fluid and magnetic field evolution 
are displayed in a conservative form in order to ease the 
comparison between analytical and numerical calculations. 
We have 

D t D + V ■ (D/3) = (2) 
AS + V • (sp -(~ + (P- B)p\ ~ + P tot l\ = (3) 
AB + V • (j3B - BP) = (4) 
D t r + V ■ (V + Aot)/3 - 09 • B)~) = (5) 

V • B = (6) 

where the conservative variables are the relativistic mass 
density D — pT, the relativistic momentum S = (T 2 ph + 
B 2 )P — (P ■ B)B, the magnetic field B and the total en- 
ergy t — r 2 ph + B 2 /4tt — Ptot respectively. The velocity 
(normalized to the speed of light c) is denoted by P while 
r~ 2 = 1 — /3 2 is the fluid Lorentz factor. The quantity 
ph stands for the enthalpy ph = pc 2 + jP/('y — 1) and 
P tot = (B 2 /V 2 + (P ■ B) 2 )/8tt + P is the total pressure. 
The derivative operator D t = dt/c. This set of equations 
is the mathematical translation of mass, momentum and 
energy conservation while the last relation ensures that the 
magnetic field is monopole-free. The characteristic speeds of 
magneto-acoustic and Alfven waves are derived by lineariz- 
ing the previous equations. In the ID configuration adopted 
hereafter the Alfven disturbance cannot propagate along the 
normal to the shock front since we assume the magnetic 
field to be parallel to the shock front. We will hence con- 
sider only magneto-acoustic disturbances even if the case of 
an Alfven disturbance is shortly discussed at the end of sec- 
tion 



2.3 (basically to show that Alfven waves are stable in 



this flow configuration). The general slow and fast magneto- 
acoustic characteristic speeds Xs,f in SRMHD are obtained 
by solving a cumbersome quartic polynomial as discussed in 



van der Hoist et al. (20081 and for which the group speed 



anisotropies are shown more clearly in Keppens & Meliani| 
(2008). In the case of a relativistic shock having a magnetic 



field perpendicular to the shock normal, the magnetic en- 
ergy density and thermal pressure are very small compared 
to the rest mass energy density. This leads to a simpler quar- 
tic polynomial which provides the expression of Xs,f 



Px and 



c 



B 1 



4-nT 2 ph 



1/2 



(') 



where the sound speed is (3 S = c,s/c = y^jP/ph and the 
Alfven speed is fi a = B y jY^fAmph. A fast magneto-acoustic 
wave propagating in the upstream medium and propagating 
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towards the shock will then display perturbed quantities as 

SB _5p = k x T 2 S/3 X _ k x r 2 sp 

p \K 



B 



+ 



B 1 



V 2 \k*\ I3 f 



(8) 



4tiT 2 p/i J 

where k x stands for the wave vector of the magneto-acoustic 
wave while /3f = Xf/c. 

To complete this study we have to account for the pres- 
ence of cosmic rays (CRs) in the shock precursor in the pre- 
vious set of equations. These supra-thermal particles carry 
electric charges that modify the Maxwell-Gauss equation 
(see also PLM09): 

V ■ E = 4tt (p p i + pgr) (9) 

where p p \ stands for the thermal plasma charge density while 
Pgr is the charge density of the CRs. The electric field is 
related to the velocity and magnetic field through Ohm's 
law stating that in ideal MHD we have 



E = -/3 x B 



(10) 



As a result the thermal plasma composing the precursor 
region is not neutral and we then have to consider the cosmic 
ray charge density in the momentum conservation equation, 
namely 



AS + V- S/3 



BB 

4tiT 2 



+ Pot I 



Pcr/3 x B 



(11) 



As the shock is perpendicular we have considered that /3 ■ 
B = 0. The electric force arises from the presence of CRs 
within the precursor and stands as an external force applied 
onto the thermal plasma. This force will then be treated as 
a source term in the numerical simulation. 

2.2 Structure of the precursor of a magnetized 
perpendicular relativistic shock wave 

The stationary Maxwell equations driving the electromag- 
netic field in the shock frame read 



E 
J 



= PzBy 



fix By 



— d x B y e z 



ppl = 



-PCR 



47T 



d x E x 



(12) 



In the steady Ampere equation above a CR current Jcr may 
be added to the right hand side. If such a current is taken 
to be parallel to the mean magnetic field, the equilibrium 
equation of motion is not modified with respect to Eq. |ll| as 
such an aligned current does not contribute to the Lorentz 
force. But CR current configurations along the shock normal 
would require a full 2D approach, which is beyond the scope 
of this work and will be considered elsewhere. 
The Faraday equation reads 



d x (p x B y ) = 



(13) 



The steady mass conservation equation links the density to 
the velocity along the normal to the shock front as 



d x (T P p x ) = 



(14) 



The stationary relativistic equation of motion in its primi- 
tive form is (see e.g. Heyvaerts & Norman (20031) 



pu ■ Vhu — ppiE H — J x B 

c 



VP 



(15) 



where u = F/3 is the spatial part of the relativistic four- 



vector velocity, such that r = (1 
remind here that h is defined as 



h = c + 



7 P 



(7 - l)c 2 



) 1/2 . Let us 



(16) 



( 7 - l)p 7 - 1 - 
The various components of the equation of motion write as 



pu x d x hu x 

-pu x d x hu y 
pu x d x hu z 



By r\ u yB y 

AtvF v ° x r 



= PCR 

= 



u z B, 



+ 



U X By 

r^O x U z 



B Hd x 



-PCR" 



(17) 

Ux By 



47rr 2 47rr 2 ^~" F 

Let us note that we added with respect to PLM09 the pos- 
sibility to have a motion u y along the mean magnetic field, 
but it is obvious that u y does not depend on x and can be 
set to in the rest of the discussion. The usual physical con- 
ditions prevailing in the interstellar medium are such that 
pc 2 > P ^ P 2 /47iT 2 h . This leads to f} a < 1, fi F < 1 and 



Under these circumstances, we can infer the balance 



of the background fluid by retaining only the zeroth order 
terms in the previous equations, namely 

2 ^ u z B y 
pc u x d x u x = — pcr _ 



pc u x d x u z = pcr 



Ux Bit 



(18) 



where we easily see that d x (u x +u z +u y ) = thus having a 
constant Lorentz factor T = r s h throughout the precursor. 
This property enables us to set a direct link between the 
mass density, the magnetic field and the velocity j3 x . Indeed, 
taking into account the previous equations, one can easily 
obtain that p/3 x and j3 x B y must remain constant along the 
normal to the shock front and hence that B y /p remains also 
constant. 

The velocity of the fluid occurring along the z direction 
is a direct consequence of the presence of charged CRs. The 
induced electromotive field is then balanced by this trans- 
verse motion of the flow. In order to assess the relative am- 
plitude of this motion, we now consider all terms appearing 
in the z-component of the equation of motion. The constant 
Lorentz factor F provides u z d x u z = —u x d x u x which leads 
to 



pc u x + 



B7, 



r 2 

1 sh 



d x i 



PCR- 



Ux By 



(19) 



471-^ r 2 h J ' r sh 

It is then easy to assess the transverse motion of the plasma 



u z (x) 



pcn(x)B y (x) 



r sh pc 2 



B 2 y(x)(rj h -i) 

A-k P {x)c*T^PI{x) 



dx (20) 



where x* is such that pcr(x > x*) = (the shock front is 
located at x — 0). The last term appearing in the denomi- 
nator is proportional to the inverse of the squared Alfvenic 
Mach number measured in the observer frame. Since only 
strong shock waves are relevant for CR acceleration (thus 
super- Alfvenic shocks), we may safely neglect this term. We 
then end up with the transverse velocity profile defined as 

By(x) 



U z (x) 



r sh p(x)c 2 



pCB,(x)dx 



(21) 
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which is dimensionless since pgr = [B/L] (L being the refer- 
ence length unit of the system). Following PLM09, we iden- 
tify this length L as the typical relaxation length of the CR 
charge density £cn = x, which leads to 



\u z \(x = 0) ~ IcrPcb.{x = 0) 



B u 



(22) 



' r shP c 2 

We define the energy density of CR in the shock restframc 
as 6cr = £cR(r s h — l)r s hp u c 2 and the charge density as 
Pcr = qecii/E t . In the previous expression, E* is the mean 
energy of the CR having their Larmor radius of the order of 
icn while q stands for the charge of a proton (E, — £cRqB y ). 
Previous relations lead to an estimate of the maximal trans- 
verse velocity, 



|«,|(z = 0) ~£c R (r sh -i) 



(23) 



where (cr is the fraction of the shock energy that is trans- 
ferred into the CR population by Fermi acceleration. The 
vertical velocity profile can now be expressed as 



vW= AM = * (F sh - 1) fi z {x) 
X{ ' &,(0) ?CR p s ^T sh fi z (0) 



(24) 



The derivative of the vertical velocity enables us to define a 
fiducial wave vector k*, namely 



fc*(a;) = 



d x u z 

r s h 



= ^CR 



(r sh - 1) pcnjx) 

^CRF sh pcr(0) 



(25) 



In the case of ultra-relativistic shock waves, we recover the 
result by PLM09, i.e. x(0) ^ £cr < 1. For mildly relativ- 
ity shocks, this assessment holds while for non-relativistic 
shocks x(0) ~ ^CRjSsh/2. The equilibrium solution used in 
this work is displayed in figure [T] 

The gradient of the magnetic field along the precursor 
induces a drift velocity/3<j and current Jd of CRs parallel to 
the shock front and perpendicular to B y and thus a force 
Jd x B y oc (B x VB) x B in the x direction. Using Eq[l3j 
Eq |25| and the conservation of the Lorentz factor of the flow 
one can deduce that |/3d| ~ £c R <C \fi z \. Hence we find the 
force produced by the CRs charge to be dominant over the 
CR drift induced force which will not be retained hereafter. 
It should be noted however that the CR gradient drift is 
a potential interesting effect to account for, once the mag- 
netic field is modified in the precursor. The drift current 
associated with the shock magnetic compression is 
difficult to infer since as already shown by Begel- 



man &: Kirk ( 1990[ ), the average number of CR shock 
front crossing is dropping for relativistic superlumi- 
nal shocks to a value close to unity. The actual am- 
plitude of the CR shock front drift current is then 
still an open issue for relativistic shocks so we dis- 
carded its effect in our calculation. The magnetic field 
amplification is investigated in the next section and the im- 
pact of the drift is postponed to a forthcoming work. 



2.3 Linear study of the streaming instability near 
perpendicular relativistic shocks 

The presence of a charge layer near the front shock will 
induce a supplementary electric force while a MHD wave 
is propagating through the layer. This force tends to am- 
plify the associated velocity disturbances that will amplify 
the magnetic disturbances through the magnetic induction 



process and so on. The system then enters an exponential 
amplification that ultimately leads to a non-linear instabil- 
ity stage. PLM09 already made a first attempt to study the 
linear stage of the instability but discarded terms associ- 
ated with the CR flow. We present in this section an alter- 
native presentation of the linear analysis of the instability 
while retaining the most significant terms leading to a larger 
growth rate of the instability. Indeed, as we show further in 
this section, PLM09 has discarded terms proportional to the 
transverse velocity f3 z arguing that its relative amplitude is 
very small. As a matter of fact, estimates achieved in (non- 
relativistic) supernovae shock waves have shown that the 
parameter £cr can reach values up to 10% (Milosavljevic & 



Nakar 2006). According to the shock equilibrium (see figure 



IJ) displayed in the previous section, this leads to a maximal 
transverse velocity of the order of c/10 in relativistic shocks 
which cannot be discarded. 

The basic assumptions made in this linear study is first to as- 
sume that both sonic (/3 S ) and Alfvenic speed (fi a ) are much 
smaller that the velocity of light. The second assumption is 
to consider perpendicular shocks, a natural shock configu- 
ration arising from the Lorentz transformation in relativis- 
tic shocks. The first assumption is induced by the physi- 
cal conditions prevailing in the interstellar medium where 
pc 2 2> -B 2 /47rFg h ^ P. The second assumption enables us 
to assess that the wave vector of any MHD wave propagat- 
ing in the precursor of a relativistic shock is mostly oriented 
along the normal to the shock front. This statement arises 
from the Lorentz transformation inducing a wavelength con- 
traction in the x direction by a factor r s h. We use then that 
k x ~S> k„ = dxUz/Ysh while in this work w6 set ky — kz — . 
In order to follow the temporal evolution of MHD waves 
entering the precursor region, we first write the set of con- 
servative RMHD equations describing the physics of the pre- 
cursor, namely 



D t (Tp) + d x {TpP x 
D t {Y 2 pc 2 fi x ) + d x {V 2 pc 2 fi 2 x 
D t {T 2 pc 2 fi z ) + d x {V 2 pc 2 fi z fi x 
DtBy + d x (B y fi x 
D t B z +d x (B z fi x 

D t {r 2 pc 2 ) + d x (r 2 pc 2 fi x 





—pcnfizBy 

PCR.fi x B y 







(26) 
(27) 
(28) 
(29) 
(30) 
(31) 



where we neglected all terms proportional either to fi a or 

Pi 



2.3.1 Magneto-acoustic wave propagation 

We consider in this study magneto-acoustic waves displaying 
magnetic perturbations along the y— direction. The induc- 
tion equation regarding B z will then be trivial as B z remains 
null during the propagation of the wave throughout the pre- 
cursor. 

The dominant force in the precursor arises from the presence 
of CR within the precursor leading to a non neutral plasma. 
The CR charge density at the shock front is provided by 



Eq.(20l 



pCR(a: = 0) 



r 2 sh pc 2 k t (x = 0) 
B v 



(32) 
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Figure 1. Initial equilibrium state of the upstream medium of an ultra-relativistic perpendicular shock (r s h = 100) with a CR electrical 
charge. In this equilibrium, the CR energy stands for £cr = 10% of the shock kinetic energy. As can be seen on the figure, the presence 
of the CR charge induces a parallel motion to the shock front (v z ) in order to balance the electrical force acting on the plasma. On 
the plots, the shock front is located at x = (outside the computational domain) and the unperturbed upstream region lies beyond 
x* = I.I^cr- Our simulations consist in injecting a magneto-acoustic wave near x„ and computing its propagation through the precursor 
from the right to the left of the box. 



Defining D x = T(D t + j3 x d x ) as a Lagrangian derivative, we 
easily get from the previous set of equations that 



D X T = 



Tpc D x j3 x 
Tpc 2 D x p z 

D X By 



By 

D x p 

p 



= pCB./3 x B y 

= -vd x p x 
= -vd x p x 



(33) 



Using a linear approach, the previous set of equations en- 
ables us to get the behavior of small perturbations as for 
instance magneto-acoustic waves propagating in the precur- 
sor. We have 



D X 5F 
D x Sp x 

D X S/3 Z 
D 5 -P 







Tk*{P Z 8By + By 



B„ 



Tk*p x 5B v 

By 

-Td x 5p x 

-rd x sp x 



(34) 



The perturbed momentum equation along the normal to the 
shock involves two terms in its rhs. Retaining only the term 
proportional to 5/3 z (thus neglecting the transverse velocity) 
leads to the instability mode found by PLM09 as combining 
the different equations provides 



r 3 kl(3 x d x s -^ 

By 



(35) 



In order to assess the linear growth rate of the magnetic 
perturbation, we use a standard harmonic description as 
D t = —ik x /3 x while d x = ik x (l — e) where e is a complex 
number whose imaginary part leads to the spatial growth 
rate of the magnetic field. The resulting growth rate of the 
instability matches the one given by PLM09, namely 



f x = Im{k x t) 



V3 
2 



k^ 



2/3 



k 



1/3 



(36) 



As mentioned earlier, the transverse velocity may be signifi- 
cant as values of £cr may reach up to 10% in the precursor. 
It seems then more realistic that \p z SB y \ 3> \B y Sp z \ regard- 
ing the propagation of magneto-acoustic waves. This then 
leads to a different equation of the evolution of the mag- 
netic perturbation, namely 



J. j-y 



-rd x D x sp x 



r 2 k,i3 z dj-^ 

By 



(37) 
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We clearly see here the source of the non-resonant instabil- 
ity which relies on the growth of the velocity perturbation 
induced by the presence of an initial magnetic perturba- 
tion. The magnetic perturbation eventually starts to grow 
thanks to magnetic induction which results into an even 
larger growth of the velocity perturbation and so on. The 
growth rate of the magnetic field will then be given by e 

2 . 



which leads in the case of small e to 



7-s = Im(k x e) 



2/3 s h 



1/2 



(38) 



(39) 



This growth rate is maximal for ultra-relativistic shocks 
(while being independent of r s h) and decreased for mildly 
relativistic and non-relativistic shocks (proportional to 
/?sh 2 in this regime). The magneto sonic waves are hence 
destabilized like in the configuration investigated by |Drury| 
& Falle (19861, but due to the electromotive force induced 



by the CR charge. 



2.3.2 Alfven waves propagation 

We have so far only considered magneto-acoustic waves and 
discarded Alfven waves since they are not prone to this insta- 
bility. Indeed, in a more general framework we can consider 
Alfven waves having wave vector k — k x e x + k y e y . The 
resulting wave is incompressible as only 8/3 z and 8B Z are 
non-vanishing contributions. The set of RMHD equations 
then reads 

rV 2 A<% = - PCR B y sp z 

r 2 pc 2 Di5l3z = pcnBySp* 

DiSB x = B y dy8p x 
DiSB y = -B y d x 8p x 

DiSB z = BydyS/3% (40) 

where Di — T(D t + fi x d x + /3 y d y ) stands for the 2D La- 
grangian derivative. We retain only terms from the pertur- 
bations 8f3 z and SB Z plus the term involved in the evolution 
of 80z, namely 8f3 x . Combining the first two equations leads 
to the evolution equation of <5/3 z 

* = (41) 

which does not have any exponentially diverging solutions. 
Thus, Alfven waves propagating through the precursor are 
found to be stable (see also PLM09) . 



3 RMHD WAVES PROPAGATING IN THE 
PRECURSOR OF A RELATIVISTIC SHOCK 

The presence of the CR electrical charge density (and the 
associated current) is the key origin of the non-resonant in- 
stability as it is able to destabilize MHD waves. We will first 
present the RMHD framework and the numerical tool we 
used to perform the non-linear simulations described in this 
work. In a second step, we present a linear description of the 



instability in a restrained but nevertheless interesting case, 
namely a transverse magnetosonic wave flowing towards the 
shock front. 



3.1 Numerical code description and initial setup 

The MPI-parallelized, Adaptive Mesh Refinement Versatile 
Advection Code (MPI-AMRVAC) is a multi-dimensional 
numerical tool devoted to solve conservative equations 



der Hoist et al. ( 2008 1; Keppens et al ( 2012 1) using finite vol- 
ume techniques and a dynamically refined grid. The MPI- 
AMRVAC package includes sets of equations as hydrody- 
namical or magnetohydrodynamical equations either in a 
classical or relativistic framework. For the simulations dis- 
played in this paper, we used a second-order Harten-Lax-van 
Leer Contact (HLLC) solver linked to a Koren slope limiter 
( Koren|1993 1 to make sure we employ a robust but accurate 
scheme. 

The AMR refinement strategy can be controlled by several 
means within the MPI-AMRVAC framework, such as by 
Richardson extrapolation to future solutions or using instan- 
taneous quantifications of the normalized second derivatives, 



or by a user controlled criterion or actually both (Keppens 



et al (2012|). For the purpose of our simulations, we simply 



choose to enforce the maximal refinement to follow the wave 
propagating throughout the precursor of the shock knowing 
its velocity as given by Eq. |7J. The base level is filled with 
blocks of equal size which can be divided into 2 D child grids 
having the same amount of grid cells than the parent grid 
(D being the dimension of the grid). The structure of the 
grid will then be similar to an octree for three dimensional 
calculations. In the context of a perpendicular relativistic 
shock, the aforementioned set of RMHD equations simpli- 
fies as B ■ (3 = 0. Assuming the magnetic field to be oriented 
along the y direction (parallel to the shock front), we obtain 
for ID simulations the following set of equations, namely 



d t D + d x (D/3 x ) = 

d t S x + d x {S x j3 x + P tot ) = -pcR0zB y 

dtSz + d x (S z /3 x ) = pcRpxBy 

d t B y + d x {B y /3 x ) = 



(42) 
(43) 
(44) 
(45) 



where we took into account the CR charge effect upon the 
momentum of matter. Considering the physical conditions 
prevailing nearby relativistic shocks, PLM09 assumed that 
both the magnetic pressure and thermal pressure can be ne- 
glected compared to mass energy density. From a numerical 
point of view this is a quite difficult feature to deal with 
as the energy equation is likely to provide numerical errors 
leading to negative thermal pressure. In order to avoid that 
problem, we decided to set a simple isothermal relation link- 
ing thermal pressure to mass density, P — pc% where c$ <C c 
is the sonic speed in the matter rest frame. The procedure to 
switch from conservative to primitive variables (p, /3,B,P) 
has then to be revisited compared to the one presented in 



van der Hoist et al. (20081 (see Appendix A). In one dimen- 



sional computations dealing with shocks having magnetic 
field parallel to the shock front, the conservation of mag- 
netic flux is naturally enforced as d x B x — 0. 
We set up the initial conditions by following the equilibrium 
presented in Sect. |2~2| ). The equilibrium is entirely controlled 
by two free parameters, namely the shock Lorentz factor r s h 
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and the amount of kinetic energy of the shock transferred 
into the CR population £cr- The parameter f3 a and p„ are 
also free but their actual value do not influence the insta- 
bility behavior provided it remains very small compared to 
unity (we specifiate their value for each simulation in the fol- 
lowing). One final ingredient is required to close the initial 
setup, namely the shape of the cosmic ray charge density. 
In order to simplify the instability study, we choose pen. to 
be constant within the precursor and vanishing outside. We 
then choose to set 

rfeH.r sh (r sh -W tanh3 /^ 

Pcr(x) = < £crB \ £ 

[0 , x x* 

(46) 

where l is a constant such that £ — 0.03^ca and B and 
p are the magnetic field and mass density at the location of 
the shock (x = 0). The computational domain ranges from 
x — where the shock front is located, to x = IAIcr- Let us 
note that we set x* — 1.1£cr to be the border between the 
precursor of the shock and the unaffected ISM medium. The 
vertical velocity /3 Z is linked to the charge density profile as 



the same technique than Meheut et al (20101, namely we 



< x ^ x, 



stated by eq.(21 1 which provides 



u z (x) 



^CR V2 V ^ 



+ In < cosh 



(47) 



The velocity along the normal to the shock can be de- 
rived from the aforementioned equilibrium, namely when 
T = r s h. The previous equilibrium is consistent with u y — 
so u x (x) = — \/r 2 h — 1 — u z {x) 2 . We also showed that the 
magnetic field components B x and B z are null. The remain- 
ing quantities (B y and p) are easily constrained by the mass 
conservation and induction equation leading to 

\l/2; 



u x (x)B y (x) 
u x {x)p(x) 



-(r 8 2 h 
-(r B 2 h 



j Pism 

i \l/2 

1) PlSM 



(48) 
(49) 



where Bism is the value of the magnetic field in the unper- 
turbed region ahead of the precursor. We choose to set the 
density of this region pism = 1 so that Pism = Pa\/i^c 2 r 2 h 
in our simulations. As an instance, we displayed in Fig.|l]) 
the variation of the various physical quantities for an ultra- 
relativistc shock having r s h = 100, £cr = 10% and P% — 



10" 



and p s 



10" 



The analytical equilibrium derived in Sect(2.2l is only an 



approximation of the actual equilibrium since we discarded 
the magnetic and thermal pressure terms involved in the 
momentum equation (Eq |15| ). The neglected pressure terms 
are very small compared to the force associated with the 
electromotive field 



By dx By 
^PCB-PzBy 



Pi 
P2» 



< 1 



d x P 
pcnPzBy 



Pi 
PL 



< 1 (50) 



In the context of the numerical equilibrium presented in 
fig.|l]), the previous ratio are equal to 10 -10 and 10 -8 respec- 
tively. From a numerical point of view this analytical equi- 
librium is not suitable for our simulations as the magneto- 
acoustic wave entering the precursor region displayed veloc- 
ity disturbances smaller than the previous ratios. We then 
have to implement a numerical approach that enforces an 
exact equilibrium to machine precision. To do so we employ 



compute the electromotive source terms in a two step fash- 
ion. During the first step, we calculate the numerical errors 
produced by the non accurate initial conditions presented in 



Sect. (2.2 I such that 



D t Si + V- (Sifr- (j| + (A • -boa) H 
+ Ptot.il) -Pcr(A x Bi) = $ t 



(51) 



where the subscript i stands for initial conditions quantities. 
Once the numerical errors are calculated, we used them into 
the second step dealing with the actual physical conditions, 
namely 



D t S + V-(s/3-(^ + (/3.B)/3)|L 
+ PtotI) = Pcr(/3 x B) + 



(52) 



Let us note that this procedure has to be repeated each time 
the adaptive mesh evolves since any modification in the grid 
hierarchy would lead to different values of <l?i. 
At the beginning of the simulation we apply the aforemen- 
tioned procedure during the first time step (and at each 
following time step where the grid hierarchy chages) and 
at the beginning of the second time step we add to the 
equilibrium physical quantities the perturbation associated 
to a magneto-acoustic wave packet in the region ahead of 
the upstream medium where no CR charge is present. More 
specifically we add perturbation to the mass density, mag- 
netic field and velocity field. The spatial extent of the wave 
packet is equal to 24ir/k x (twelve wavelengths) as can be 
seen on the upper left plot of Fig.|2|. The perturbation is 
coherent with a magneto-acoustic wave whose wave vector 
is k = k x e x in the shock rest frame. The detailed expression 
of the implemented perturbation is provided by Eq.|8|. The 
boundary conditions are simply continuous boundary where 
all physical quantities are copied in the ghost cells. 



3.2 A reference simulation of the instability 

3.2.1 The basic stages of the instability 

The initial perturbation entering the CR charge dominated 
region is prone to the additional force arising from the CR 
charge density. As it can be seen from Eq.Q, the initial 
perturbed velocity 8p x is very small compared to 5B y /B y 
so the growth of the velocity perturbation will not induce 
any amplification of the magnetic disturbance right from the 
start. During this preliminary stage, it is possible to infer 
the growth of the velocity perturbation by looking at the 
linearized SRMHD equations. Taking into account the mass 
density conservation and the induction equations leads to 



d t 5p ~ -p x d x 5p - d x 6px 
d t 8B y ~ ~p x d x 5B y — d x 8f3 x 



(53) 
(54) 



which are fully consistent with a pure advection motion as 
long as Sp x and SB y /B y remain synchronized. On the other 
hand, the linearized equation of motion lead to the following 
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Figure 2. Snapshots of the reference simulation presented in Sect. |3~2.2| . We have displayed the parallel velocity perturbation together 
with the magnetic perturbation at locations related to the various stage of the instability. At first (top left panel) the initial magneto- 
acoustic perturbation is entering the cosmic ray dominated region then (top right) the velocity perturbation is increasing while the 
magnetic perturbation remains unchanged. The third step (middle left) shows the beginning of the magnetic field amplification when 
both perturbations start to be out of phase. The fourth stage (middle right) shows that the magnetic amplification starts to modify 
the shape of the magnetic perturbation because of the phase shift (being different from tt/2). The fifth stage (bottom left) shows that 
once the magnetic perturbation is modified, it then affects the shape of the velocity perturbation. The final stage corresponds to a fully 
non-linear stage where the magnetic field is completely dominated by the magnetic disturbance. The figures also display a time sequence 
evolving from the upper left (perturbation entering the precursor at X max /£cR = 1-1) to lower right (perturbation reaching the shock 
front at X = 0) with timescales given by (X max — X)/f} s f l , where X max / ^cr = 1.1- 



relations The presence of the cosmic ray charge density alters the 

nature of these equations since they are no longer consistent 
with a pure advection motion. The amplitude of the velocity 
a r R fl to disturbances increases while the wave is propagating and 

Otopx + PxOopx + j^~o x -^— ~ — ?cr^ — p- (55) we can see that the velocity perturbation is then (assuming 

the disturbances propagates at speed very close to /3 X which 



r s \ By - ^£cnB y 
/3 F a 5B y ^ f p x 8B y 



d t 5/3 z + /3 x d5f3 z + ^4j-d x — - ~ £cr (56) remains almost constant since /3 Z is very small at the outer 

T sh By ecnBy 
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Figure 3. The mass density and magnetic perturbations at the 
final stage of the simulation. The system has encountered the 
saturation limit expected from the amplification of a magneto- 
acoustic wave, namely the regime where Sp/p — > —1. We can see 
that the initial wave has been amplified to a point where the 
density of the medium exhibits very sharp variations that may be 
associated with shock fronts within the precursor. 



edge of the CR charge dominated region) 
Pl{x) 8 By 



8/3 z {x) 



8/3 x (x,) + 

By 



2/?sh Bi 



'-{1+1%) 



(57) 



where we used the fact that prjR. — (r sri pc 2 d x u z ) / B y (see 



Eq 18 1. During this first stage, only the velocity disturbances 
are amplified while the magnetic perturbation remains un- 
changed while propagating. Nevertheless the growth of the 
velocity disturbance has an impact on the advection veloc- 
ity of the magnetic perturbation as can be seen from the 
linearized induction equation, namely 



dt5B y + [/3 X + 



r 2 

1 sh 



+ 



(1 + /3 2 ) )d x 6B v =0 (58) 



The advection velocity of the magnetic disturbance is 
clearly depending on the transverse velocity of the flow in- 
duced by the CR charge density. Since the advection speed 
of the velocity disturbances remains constant, the magnetic 
and mass density perturbations begin to be out of phase with 
the velocity disturbances which leads to an amplification of 
both magnetic and mass density perturbations. The wave 
is then entering a second stage where the instability arises 
from both magnetic and velocity exponential growth. In or- 
der to assess the location of the transition between the two 
stages, we may look at the previous equation where falling 
back to the general expression, we have 



d t 8By + p x d x 8B y + B y d x 8p x 



d8B. 



y -+B v d x 5p x =0 (59) 



The transition between the pure advection stage and the 
magnetic field amplification stage occurs when the veloc- 
ity perturbation is sufficiently large to generate a variation 
of the magnetic perturbation of the order of SB y /B y when 
crossing the entire CR charge dominated region, namely 
when A5By/B y ~ SB y /B y . This leads to a velocity per- 



Figure 4. The maximum amplitude of the magnetic perturbation 
as a function of the location of the perturbation. The exponential 
growth of the perturbation matches the linear estimate of the 
instability. 



turbation threshold 



toTians /3shA 8 By ,0 Trans 

S/3 X ~ 5-^ and /3 Z 

«CR By 



^CR 



(60) 



where A is the typical wavelength of the magneto-acoustic 
wave. It is obvious that shorter wavelength waves will ignite 
the magnetic field amplification earlier as the threshold re- 
lation is met earlier in the wave propagation. 
At a certain point, a third stage is expected as the phase 
shift between magnetic and velocity disturbances continues 
to increase during the second stage. Once the offset is suf- 
ficiently large, the magnetic perturbation will modify the 
shape of the velocity perturbation since its maximal growth 
will occur in parts of the wave where the velocity disturbance 
is not maximal. The modification of the velocity perturba- 
tion will then alter the magnetic disturbance profile and so 
on. This third stage is consistent with a fully non-linear dy- 
namics. 



3.2.2 A reference numerical simulation in an 
ultra-relativistic shock 

In order to illustrate the aforementioned stages of the 
instability, we designed a reference simulation matching all 
physical conditions assumed in the linear analysis, namely 
where B 2 /4vr ~ 10 2 P ~ 10~ 8 pc 2 . We set up an ultra- 
relativistic shock with r sn = 100, ensuring that the initial 
magneto-acoustic perturbation is such that 8f3 x <C SB y /B y . 
In order to get a large growth rate, we set the ratio of 
CR. energy to shock energy <fcR 



10% ( Milosavljevic fc 



Nakar 2006). The various physical quantities are initially 
set according to the equilibrium presented in Sect. |2~2| ) and 
displayed in Fig.|l]). The numerical equilibrium is achieved 
within machine precision according to the procedure 
presented in the previous section. We start the simulation 
with the previous setting and add a magneto-acoustic wave 
according to Eq.|8| at the outer right hand side border of 
the upstream region, namely where pen = 0. The wave 
vector of the perturbation is such that k x £cR — 2.5 x 10 4 . 
In order to capture the dynamics of the instability, we used 
the adaptive mesh refinement technique. The base level of 
the simulations has 1700 cells while 10 sub-grid levels were 
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used, reaching an effective resolution of near 8 x 10- 7 ^cr. 
Computing the crossing from right to left of the entire 
region by the wave is quite computationally challenging 
since it requires several millions iterations. 

We displayed in Fig.|2| various snapshots of both the 
velocity and magnetic perturbations from the numerical cal- 
culation. Any perturbed quantity A is computed as SA — 
A(t) — A(t = 0) where A(t = 0) is the quantity as defined by 
the initial equilibrium. In the top left panel we displayed the 
initial perturbation setup where we can see the ratio of the 
relative amplitudes is 10 s = T^/Pf. On the top right panel, 
the wave has entered the charge layer but is located in the 
region where pcr(^) is increasing. The magnetic perturba- 
tion remains unchanged while the velocity perturbation has 
grown already. The amplitude of the velocity perturbation 



matches the assessment from Eq.(57l very well. The mid- 



dle left plot shows the wave at a stage where the magnetic 
perturbation is starting to grow. At this point, the mag- 
netic perturbation is slightly out of phase with the veloc- 
ity component as expected. Let us note that the magnetic 
field amplification occurs when the ratio between the ve- 
locity perturbation to the magnetic perturbation is of the 
order of 2 x 10~ 4 which is close to the set wavelength of the 
magneto-acoustic wave, namely A = 2.5 x 10 _4 ^cr. This 
is in agreement with the magnetic amplification threshold 



identified in Eq. ( 60 1 



The middle right panel shows the two perturbations 
at a later stage where the phase shift between the two 
components becomes quite significant. This shift induces 
a modification of the perturbation profile. The bottom 
left plot shows the status of the system while entering 
the non-linear regime, namely when SB y /B y is no longer 
<C 1. The perturbation profiles are completely distorted by 
non-linear effects. Finally the bottom right plot shows the 
system at its final stage, namely where saturation occurs. In 
this final snapshot the magnetic perturbation has reached 
an amplitude larger than the original magnetic field, leading 
to a modest but significant magnetic amplification (up to 
5 times the initial magnetic field). As already mentioned, 
mass density and magnetic field are controlled by similar 
equations so it is expected that the mass density reaches a 
similar state than the magnetic field. As shown in Fig.(|3|, 
it is verified by the numerical simulation since both mass 
density and magnetic field perturbations are very similar. 
The magnetic field fluctuations significantly exceed those 
of the initial perturbation such that SB/B > 0, meaning 
that the magnetic flux has increased and that the magnetic 
field has been amplified. The amplification of the magnetic 
field and mass density leads to regions where Sp/p — > — 1, 
namely regions whose mass density tends toward zero. In 
our simulations, this regime ultimately leads to a numerical 
instability, since true zero densities cannot be handled by 
the MHD code, and fiducial fixes to avoid negative densities 
become activated. 

The temporal evolution of the magnetic perturbation 
amplitude is displayed in Fig.Q. The spatial growth of the 
magnetic disturbance does not show a constant exponential 
growth during the precursor crossing. This was expected as 
the spatial growth rate inferred from linear analysis involves 
functions x( x ) an d k* (x) depending of the distance x from 



the shock front. We displayed on Fig.Q the growth rate 
predicted from Eq.(39l. Both the numerical and theoretical 



growth rates match quite well during the first part of the 
simulation while the perturbation remains very small com- 
pared to the mean magnetic field. Once the perturbation has 
grown, the numerical growth rate becomes larger than the 
linear prediction. This is not surprising since at this stage 
SBy/By has reached a relative amplitude of a few percents 
and the linear assumption no longer holds. We can observe 
at the end of the simulation that the saturation process oc- 
curs leading to a slight flattening of the growth rate. 



3.2.3 Non-resonant streaming instability; from 
ultra-relativistic to non-relativistic shocks 

We have run similar simulations to the one presented in the 
previous section but with various shock velocity and Lorentz 
factors. The growth rate predicted by the linear analysis is 



Jx(x) = Im{k x e) 



k*{x)x{x)k x 
2/Ssh 



1/2 



(61) 



where x( x ) an< A k*(x) are functions such that the spatial 
growth reads 



7*(z = 0) 



CcR(F sh - 1) 



2/3 s h^CR 



1/2 



(62) 



In the case of ultra-relativistic shocks, we easily get that the 
growth rate is independent of the shock velocity as 



7a (a: = 0)~£ 



CR 



f k x \ 
\2£crJ 



1/2 



(63) 



The mildly relativistic shock regime is more cumbersome. 
Since we can assume /3 s h to remain close to unity in this 
regime, we see that the growth rate is a slowly decreasing 
function of . 

For non-relativistic shocks, the spatial growth rate reads 



lx (x = 0) ~ £ c 



3/2 



1/2 



(64) 



which is dependent on the shock velocity. According to the 
previous expression, we thus expect the instability efficiency 
to decrease as the shock velocity decreases. 

We have displayed in Fig.|5| the spatial growth rates 
obtained from various simulations with different Lorentz fac- 
tors. The top curve (displayed with circles) corresponds to 
two simulations with r s h = 100 and 50: both curves are 
quasi identical proving that for ultra-relativistic shocks, the 
spatial growth rate is independent of the shock velocity or 
Lorentz factor (the linear growth rate relative difference be- 
tween the two simulaton is less than 0.7% so that we only 
plotted the linear growth rate related to the r s h = 100 sim- 
ulation). A similar simulation but with T s h — 10 is plotted 
using stars. It exhibits a similar behavior than the previous 
curves but with a slightly smaller slope, as expected since 
this simulation is consistent with a mildly relativistic shock. 
The same behavior is observed for shocks having F s h = 5, 3.5 
and 2. Let us note that all curves match well the linear es- 



timate from Eq.(39l as long that the relative magnetic per- 



turbation remains small compared to unity. The last plot 
(displayed with squares) stands for a simulation of a shock 
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Figure 5. Same plot than in Fig.|4]| but for various Lorentz factors ranging from = 100 down to r s h = 1.5. The spatial growth rates 
corresponds to similar simulations than the one presented in Sect, |3~2,2| but for various shock velocity (F s ^ = 100, 50, 10, 5, 3.5, 2, 1.5). 
The linear growth rate provided by Eq.l |39| is displayed using dashed lines. We clearly see that once entering the non-relativistic regime, 
the growth rate of the magneto-acoustic waves tends to zero thus no instability is observed. The numerical curves differ from the linear 
growth ratesat the beginning of the propagation since the instability is not at work right from the start of the simulations but instead 
has to get through the various stages described by Sect.l pO} . We then adjust the linear curves to fit the numerical calculations by tuning 
the linear original amplitude of the instability. 



whose velocity is /3 s h — 0.74 (T s h = 1.5). The slope of the 
plot is stalling in that case and almost no wave amplifica- 
tion occurs during the crossing of the precursor. We then 
confirm what was expected from the linear analysis, namely 
that this instability is relevant only for ultra-relativistic and 
even mildly relativistic shock precursors. 

3.2.4 Investigating the wave number dependence 

We have performed several simulations of upstream media 
of ultra-relativistic shocks having r s h = 100 for various 
magneto-acoustic wavenumbers ranging from A = 6.25 x 
W~ 5 £cr to A = 10 _3 ^cr- The growth rates resulting from 
these simulations are displayed in Fig.(|6|. We once again 
see that the numerical results are consistent with the linear 
growth rate since the growth rate of the perturbation is in- 
creasing with the wave vector of the wave k x . More precisely, 
1/2 

we see that the k x dependence of the linear prediction is 
confirmed as the linear growth rates (displayed as dotted 
lines) fit the numerical growth rates well as long as the rel- 
ative perturbation is small compared to unity. This numeri- 
cal experiment confirms the theoretical prediction that these 
small scale waves are likely important contributors to mag- 
netic turbulence in the shock precursor. 
In order to evaluate the maximal growth rate expected from 



such instability, we have considered waves with the small- 
est wavelength possible, namely waves having k x ~ £mhd 
where ^mhd is the scale limit for the validity of MHD. The 
definition of ^mhd is PaRl, p i where /3 a = V a /c and Rl, p i 
is the Larmor radius of the thermal protons. On the other 
hand the size of the precursor £cn is of the order of the Lar- 
mor radius Rl* of the most energetic CR generated by the 
shocks. The largest growth rate in an ultra-relativistic shock 
will then be 

( R \ 1/2 ( Y \ 1/2 

7 r^cH \2R^i) ~ ic A^) (65) 

where F, is the Lorentz factor of the most energetic CRs 
generated by the shock. For standard values of density and 
magnetic field occurring in the ISM, one gets /3 a ~ 10 -4 . 
The previous growth rate is much larger than unity as long 

as 

Ccr > I — p J (66) 

which for an ultra-relativistic shock Y s h = 100 and CRs with 
maximal energy E, = 10 15 EeV is equal to £cr ^ 10~ 4 . 
It seems then most likely that efficient wave amplification 
may be achieved in the precursor of the shock. However it 
is too early to raise definite conclusions here, as only multi- 
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Figure 6. Same plot than in Fig.Q but for various wavenumber ranging from k x £cn = 6.25 X 10 3 up to k x £cR = 10 5 . The spatial 
growth rates correspond to similar simulations than the one presented in Sect. ( [3T2.2[ | except for the wavenumber of the perturbation 
(k x £ C n = 6-25 X 10 3 , 1.25 X 10 4 , 2.5 X 10 4 , 5 X 10 4 and 10 5 ). The linear growth rates are obtained from Eq.f39l. The numerical growth 



rates obtained from simulations follow the predicted linear growth rate which is proportional to k. 



1/2 



dimensional calculations can catch the diversity of the source 
of free energy and provide a hint of the dominant MHD 
instability. We insist also on the issue that our work can- 
not provide any answer to the acceleration of high energy 
CRs in relativistic shocks yet, as again 3D calculations are 
mandatory for that purpose. The fact that the instability 
investigated here produces perturbations at scales A < tl 
makes it necessary to also produce longer wavelength per- 
turbations. Similar issues have been discussed in the context 
of supernova remnant shocks ( Schure et al|2012 l. 



4 SUMMARY AND CONCLUDING REMARKS 

In this paper we have provided an analysis of the non- 
resonant streaming instability occurring in the precursor 
of astrophysical shock waves. Our investigation consists of 
a linear analysis complemented by numerical relativistic 
MHD simulations. Two assumptions have been made to 
perform our work, namely that the MHD waves have 
a wavevector perpendicular to the shock front and a 
magnetic field to be parallel to the shock front. This 
configuration is likely to occur near ultra-relativistic 
shock fronts because of the properties of the Lorentz 
transformation. The non-resonant streaming instability is 
based on the presence of CRs in the precursor that lead 
to an electrically charged thermal plasma. The resulting 
electromotive force on the thermal plasma can destabilize 



MHD waves propagating through the precursor of the shock. 

The outcome of our linear analysis is presented in 
Sect.|2| and shows that under the aforementioned assump- 
tions, magneto-acoustic waves will be destabilized by the 
electromotive force and can reach a saturation regime if 
its wave vector k x is large enough compared to the inverse 
of the size of the precursor £cr- The linear calculation 

1/2 

predicts a k x dependence regarding the growth rate of the 
incoming wave. It also predicts that the shock velocity can 
influence the actual value of the growth rate as the growth 
rate is also proportional to (r s h — I) /T^fi^ . This factor 
has a maximal value of unity for ultra-relativistic shocks 
while it is decreasing with shock velocity ,9 s h as one enters 
the non-relativistic shock regime. Ultra-relativistic shocks 
are then the best candidates for efficient non-resonant 
streaming instability. 

We have performed numerical relativistic MHD 
simulations in order to follow the propagation of a 
magneto-acoustic wave through the precursor of a magne- 
tized shock. Under the same framework than in the linear 
analysis, we were able to depict the whole propagation of 
the wave from its entrance in the precursor to the non linear 
saturation regime. We have illustrated the basic stages 
of the instability and were able to compare the growth 
rate provided by numerical simulations with the computed 
linear growth rate. The numerical simulations match the 
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growth rate predicted by linear theory and confirm that 
under our assumptions, non-resonant streaming instability 
is more efficient in the context of ultra-relativistic shocks 
than in the non-relativistic case. The growth rate of the 
instability is found to be proportional to k\ {k x being the 
wavenumber of the perturbation). Numerical simulations 
also highlight the non-linear growth of the instability in 
the final stages of the propagation of the wave. During 
that stage, the shape of the initial perturbation has been 
completely distorted while the amplitude of the magnetic 
perturbation has increased with 8B/B of the order of a few 
units. This stage corresponds to the saturation regime of 
the instability since the density perturbation is following 
the same behavior than the magnetic perturbation. Ul- 
timately though, numerical errors occur when in the full 
saturation regime, we get true local vacua with Sp/p — > — 1. 
We found that the precursor of the shock can easily be 
populated by a train of shocks produced by the density and 
magnetic spikes. These shocks can actively participate to 
the local heating process of the background plasma as well 
as possibly to further particle acceleration and transport. 
This aspect will be investigated in a future work. We also 
confirmed that Alfven disturbances are stable in this ID 
flow configuration. 

We restrict ourselves in this work to one-dimensional 
simulations and linear calculations of magneto-acoustic 
waves amplification. This choice was motivated by the nat- 
ural configuration expected from an ultra-relativistic shock 
where the Lorentz transform leads to perpendicular mag- 
netic fields and MHD waves having k x k y ,k z . Compared 
to PLM09, we identified a new efficient instability regime 
that we believe is the most dominant mode among all in- 
stability regimes triggered by the presence of CR electrical 
charge density. Only considering multi-dimensional calcula- 
tions will permit to isolate the dominant instability that can 
provide magnetic amplification. Also 3D simulations are nec- 
essary to account properly for the particle transport and de- 
rive the maximum CR energies in relativistic shocks. Hence 
we make it clear that the question of the dominant MHD in- 
stability and maximum CR energy in relativistic and mildly 
relativistic shocks can not be discussed within this work. A 
further step will be to deal with two-dimensional simulations 
of the system but this time taking into account the effect of 
the electrical current driven by the CR flow. This external 
current may be a promising element to the shock structure 
since it may have the ability to amplify Alfven waves and 
thus provide higher SB/B ratios. 
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APPENDIX A: ISOTHERMAL SRMHD: 
SWITCHING FROM CONSERVATIVE TO 
PRIMITIVE VARIABLES 

The calculations presented in this paper involves physical 
conditions where a large contrast between physical quanti- 
ties occur. Indeed, in the context of ultra-relativistic MHD 
shocks, thermal pressure and magnetic pressure existing in 
the ISM are sensed as very small in the shock frame because 
of the Lorentz transformation. For instance, the simulation 



presented in Sect.(3.2| exhibits thermal pressure 10 times 



smaller than the kinetic energy density of matter. Dealing 
with such contrast is tricky and requires a specific approach. 
In order to avoid having negative thermal pressure in our 
simulation we add an extra relation linking thermal pressure 
to density such that P — pc% where eg ^ is the con- 
stant sound speed in the shock frame. Indeed, performing 
SRMHD simulations requires very frequent switches from 
conservative (D, S ,r, B) to primitive variables (p,v, P, B). 
In order to achieve this transformation, it is convenient 
to use auxiliary variables such as the Lorentz fac tor V or 
relativistic enthalpy £ = F 2 (pc 2 + 7-P/(7 — 1)) (see van der 



Hoist et al. (2008)) 



The major issue in switching variables comes from the 
fact that computing auxiliary variables leads to complex 
equations. In the case of isothermal SRMHD, the two aux- 
iliary variables read as 

l_ (S + r 1 (S-B)B) 2 

r 2 



£ = rz? l 



(C + s 2 ) 

7 4 
7 — 1 c 2 



= VDF 



(Al) 
(A2) 



where B 2 denotes the magnetic pressure (magnetic perme- 
ability of vacuum set to unity) and J- is a constant. Solving 
this systems leads to a fourth order polynomial in V 

F 4 D 2 T 2 + 2r 3 DTB 2 + T 2 (B 4 - S 2 - D 2 T 2 ) - (A3) 



21' j B 2 DT+ {S D ^ Y 



4 (S ■ ByB 2 

B - D 2 T 2 - ° 



Hopefully, a single root of this polynomial corresponds to the 
actual value of the Lorentz factor. This is easy to show in the 
present framework as we can assume B 2 <C D and T a 1. 
The momentum then becomes S 2 ~ T 2 D 2 (T a being the 
root of the polynomial) and the polynomial can be written 
as 



f(T) = F 4 D 2 + 2T 3 DB 2 



2T B D + 



D 



F 2 (F 2 + 1)D 2 - 

r 4 D 2 -r 2 (r 2 + i)D 2 (A4) 



The relevant values of the Lorentz factor are such that T ^ 1 
so we have to prove that the polynomial only has one root 
in the range [l,+oo[. In that context, It is obvious that 
/(r = 1) ~ —D 2 < while 



lim f(F) = r 4 L> 2 > 



The derivative of f(T) is 
dT 



1) 



(A5) 



(A6) 



which is negative if T < r cri t = \J{? 2 + l)/2 < F a and 
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positive otherwise. The function / is then decreasing and 
remaining negative for F € [l,F cr it] and strictly increas- 
ing for r 6 ]r c rit,+oo[. It is then obvious that only one 
root of the polynomial lies in the interval [1, +oo[ as long as 
B 2 ,P <C D. Finding the root of function /(F) is achieved 
using Newton-Raphson algorithm displaying a quartic effi- 
ciency (see van der Hoist et al. (20081). 
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